/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */

//
// 文件: GPUConcreteKernels.h
// 创建时间: 2025年10月4日
// 说明: GPU钢筋混凝土本构核函数声明
// 作者: OpenSees GPU加速项目组
//
// 功能:
//   - 定义GPU混凝土材料参数和状态结构
//   - 声明塑性损伤混凝土本构算法（Lee & Fenves 1998）
//   - 支持拉伸/压缩损伤、塑性应变演化
//

#ifndef GPUConcreteKernels_h
#define GPUConcreteKernels_h

#include <cuda_runtime.h>

//==============================================================================
// GPU混凝土材料参数结构
//==============================================================================

/**
 * @brief 塑性损伤混凝土材料参数（PlasticDamageConcrete3d）
 *
 * 基于Lee & Fenves (1998)塑性损伤模型
 * 参考: Plastic-Damage Model for Cyclic Loading of Concrete Structures
 */
struct PlasticDamageConcreteParams {
    // 弹性参数
    double E;              // 弹性模量 (Pa)
    double nu;             // 泊松比

    // 强度参数
    double ft;             // 抗拉强度 (Pa)
    double fc;             // 抗压强度 (Pa)

    // 塑性参数
    double beta;           // 塑性流动参数 (0~1)

    // 损伤演化参数
    double Ap;             // 拉伸损伤参数
    double An;             // 压缩损伤参数
    double Bn;             // 压缩损伤指数
};

/**
 * @brief FAReinforcedConcrete材料参数
 */
struct FAReinforcedConcreteParams {
    double E;              // 弹性模量
    double nu;             // 泊松比
    double fc;             // 抗压强度
    double ft;             // 抗拉强度
    double epsc0;          // 峰值压应变
    double epscu;          // 极限压应变
};

/**
 * @brief ConcreteL01材料参数
 */
struct ConcreteL01Params {
    double E;              // 弹性模量
    double fc;             // 抗压强度
    double epsc0;          // 峰值应变
};

/**
 * @brief GPU混凝土材料状态变量
 * 存储每个积分点的历史变量
 */
struct ConcreteState {
    // 损伤变量
    double dp;             // 拉伸损伤变量 (0~1)
    double dn;             // 压缩损伤变量 (0~1)

    // 损伤阈值（历史最大值）
    double rp;             // 拉伸损伤阈值
    double rn;             // 压缩损伤阈值

    // 塑性应变（6个分量：εxx, εyy, εzz, γxy, γyz, γzx）
    double eps_p[6];

    // 有效应力（6个分量）
    double sige[6];

    // 辅助变量
    double stateVars[8];   // 预留的状态变量
};

//==============================================================================
// GPU混凝土本构算法声明
//==============================================================================

/**
 * @brief GPU塑性损伤混凝土应力更新核函数
 *
 * 算法流程（Lee & Fenves 1998）：
 * 1. 塑性部分：
 *    - 计算试探有效应力 σ̃_trial = σ̃_prev + Ce : Δε
 *    - 应力分解为正负部分（代数法）
 *    - 检查屈服条件（负应力不变量）
 *    - Return Mapping塑性修正
 *
 * 2. 损伤部分：
 *    - 计算等效应力 τ+ (拉伸), τ- (压缩)
 *    - 更新损伤阈值 rp, rn
 *    - 计算损伤变量 dp, dn
 *    - 最终应力 σ = (1-dp)σ+ + (1-dn)σ-
 *
 * 3. 切线刚度：
 *    - 计算一致性切线矩阵 Cep
 *
 * @param strain_increment 应变增量 Δε [6]
 * @param strain_total 总应变 ε [6]
 * @param params 材料参数
 * @param state 材料状态（输入输出）
 * @param stress 输出应力 σ [6]
 * @param D_tangent 输出切线刚度矩阵 [36] (6x6按行存储)
 */
__device__ void computePlasticDamageConcrete3D(
    const double* strain_increment,
    const double* strain_total,
    const PlasticDamageConcreteParams& params,
    ConcreteState& state,
    double* stress,
    double* D_tangent
);

/**
 * @brief FAReinforcedConcrete应力更新
 */
__device__ void computeFAReinforcedConcrete(
    const double* strain_increment,
    const FAReinforcedConcreteParams& params,
    ConcreteState& state,
    double* stress,
    double* D_tangent
);

/**
 * @brief ConcreteL01应力更新
 */
__device__ void computeConcreteL01(
    const double* strain_increment,
    const ConcreteL01Params& params,
    ConcreteState& state,
    double* stress,
    double* D_tangent
);

//==============================================================================
// GPU辅助设备函数声明
//==============================================================================

/**
 * @brief 计算弹性刚度矩阵（3D各向同性）
 *
 * Ce = K * I ⊗ I + 2G * Id
 * 其中：K = E/(3(1-2ν)), G = E/(2(1+ν))
 *
 * @param E 弹性模量
 * @param nu 泊松比
 * @param Ce 输出：弹性刚度矩阵[36] (6x6)
 */
__device__ void computeElasticStiffness3D(
    double E,
    double nu,
    double* Ce
);

/**
 * @brief 应力分解函数（代数法）
 *
 * 将应力张量分解为正负部分：
 * σ+ = (σ + |σ|)/2
 * σ- = σ - σ+
 *
 * @param stress 输入应力[6]
 * @param stress_pos 输出正应力[6]
 * @param stress_neg 输出负应力[6]
 * @param Q_pos 输出投影矩阵[36] (6x6)
 * @param Q_neg 输出投影矩阵[36] (6x6)
 */
__device__ void stressDecomposition(
    const double* stress,
    double* stress_pos,
    double* stress_neg,
    double* Q_pos,
    double* Q_neg
);

/**
 * @brief 计算应力不变量（八面体应力）
 *
 * σ_oct = (σxx + σyy + σzz)/3
 * τ_oct = sqrt(2J2/3)
 *
 * @param stress 应力[6]
 * @param sigoct 输出：八面体正应力
 * @param tauoct 输出：八面体剪应力
 */
__device__ void computeStressInvariants(
    const double* stress,
    double& sigoct,
    double& tauoct
);

/**
 * @brief 计算应力范数
 *
 * ||σ|| = sqrt(σ11² + σ22² + σ33² + 2σ12² + 2σ23² + 2σ13²)
 *
 * @param stress 应力[6]
 * @return 应力范数
 */
__device__ double computeStressNorm(const double* stress);

/**
 * @brief 矩阵向量乘法（6x6矩阵 × 6向量）
 *
 * result = matrix * vector
 *
 * @param matrix 矩阵[36] (6x6按行存储)
 * @param vector 向量[6]
 * @param result 输出向量[6]
 */
__device__ void matrixVectorMult6x6(
    const double* matrix,
    const double* vector,
    double* result
);

/**
 * @brief 向量点积（6维）
 *
 * @param v1 向量1[6]
 * @param v2 向量2[6]
 * @return 点积结果
 */
__device__ double vectorDot6(const double* v1, const double* v2);

/**
 * @brief 矩阵求逆（6x6对称矩阵）
 *
 * 使用Gauss-Jordan消元法
 *
 * @param matrix 输入/输出矩阵[36]
 * @return 0-成功，-1-奇异矩阵
 */
__device__ int matrixInverse6x6(double* matrix);

/**
 * @brief 张量积运算（向量×向量 → 矩阵）
 *
 * result = v1 ⊗ v2
 *
 * @param v1 向量1[6]
 * @param v2 向量2[6]
 * @param result 输出矩阵[36] (6x6)
 */
__device__ void tensorProduct6x6(
    const double* v1,
    const double* v2,
    double* result
);

//==============================================================================
// C接口声明（用于从C++代码调用）
//==============================================================================

#ifdef __cplusplus
extern "C" {
#endif

/**
 * @brief 初始化材料状态变量
 *
 * 根据材料参数初始化ConcreteState结构体，特别是损伤阈值：
 * - rp = ft / E (拉伸损伤阈值)
 * - rn = fc / E (压缩损伤阈值)
 *
 * @param numGaussPoints 积分点数量
 * @param d_elementData GPU单元数据数组（包含材料参数）
 * @param d_states GPU材料状态数组（输出）
 * @param numElements 单元数量
 * @param gaussPointsPerElement 每个单元的积分点数量
 * @param stream CUDA流
 * @return 0-成功，-1-失败
 */
int initializeMaterialStates(
    int numGaussPoints,
    const void* d_elementData,
    ConcreteState* d_states,
    int numElements,
    int gaussPointsPerElement,
    cudaStream_t stream
);

/**
 * @brief 启动混凝土材料kernel（测试接口）
 *
 * @param numGaussPoints 积分点数量
 * @param d_strainIncr GPU应变增量数组
 * @param d_strainTotal GPU总应变数组
 * @param d_params GPU材料参数数组
 * @param d_states GPU材料状态数组
 * @param d_stress GPU输出应力数组
 * @param d_tangent GPU输出切线刚度数组
 * @param stream CUDA流
 * @return 0-成功，-1-失败
 */
int launchConcreteKernel(
    int numGaussPoints,
    const double* d_strainIncr,
    const double* d_strainTotal,
    const PlasticDamageConcreteParams* d_params,
    ConcreteState* d_states,
    double* d_stress,
    double* d_tangent,
    cudaStream_t stream
);

#ifdef __cplusplus
}
#endif

// ===================================================================
// 设备函数实现（inline实现以支持跨编译单元调用）
// ===================================================================
#ifdef __CUDACC__
#include "GPUConcreteKernels_impl.cuh"
#endif

#endif // GPUConcreteKernels_h
